﻿using com.foxmail.wyyuan1991.Assignment.Solver.Data;
using System;
using System.Collections.Generic;

namespace com.foxmail.wyyuan1991.Assignment.Solver
{
    public class PenaltySolver
    {
        #region 属性
        /// <summary>
        /// 收敛阈值
        /// </summary>
        public double Epsilon { get; set; }//epsilon = 0.0000001m;
        /// <summary>
        /// 罚因子
        /// </summary>
        public double Sigema { get; set; }//Sigema = 99m;
        /// <summary>
        /// 罚值/拉格朗日乘子 调整阈值
        /// </summary>
        public double Beta { get; set; }//beta = 0.001m;
        /// <summary>
        /// 默认步长
        /// </summary>
        public double Alpha { get; set; }//alpha = 2m;
        /// <summary>
        /// 迭代上限
        /// </summary>
        public int IterationUpbound { get; set; }
        /// <summary>
        /// 最短路搜索算法
        /// </summary>
        public IRouteSearchAlgorithm SA { get; set; }
        /// <summary>
        /// 线性搜索算法
        /// </summary>
        public ILineSearchSolver LS { get; set; }
        /// <summary>
        /// 全有全无分配算法
        /// </summary>
        public AoNSolver AS { get; set; }
        #endregion

        #region Events
        //一次循环完成事件
        public event ALSIterationFinishedEventHandler AFEvent;
        private void RasieALSIterationFinishedEvent(int i, double p, double e)
        {
            if (AFEvent != null)
            {
                AFEvent(this, new ALSIterationFinishedEventArgs(i, p, e));
            }
        }
        #endregion

        /// <summary>
        /// 罚函数法
        /// </summary>
        /// <param name="net"></param>
        /// <param name="f"></param>
        /// <returns></returns>
        public IFlowPattern Solve(INetwork net, IDemand f, List<SideConstraint> SCList)
        {
            FWSolver s = new FWSolver()
            {
                kesi = 0.000001,//设置收敛阈值
                SA = new DijstraRouteSearcher() { Network = net },
                AS = AS,
                LS = new DichotomyLineSearcher() { Kesi = 0.00000001 }
            };

            //初始化拉格朗日乘子
            //SCLagrangeMultiplier lm_s = new SCLagrangeMultiplier();//这里的拉格朗日乘子是一维的，因为边约束的形式是link capacity.
            //SCLagrangeMultiplier lm_u = new SCLagrangeMultiplier();
            
            //初始假设 对偶变量为0
            FlowPattern fd;
            FlowPattern fd2 =  s.Solve(net, f) as FlowPattern;
            int iter = 0;
            //③罚值是否超限
            while (iter < IterationUpbound)
            {
                if (iter++ > IterationUpbound) break;
                fd = fd2; fd2 = null;
                //②内嵌F-W
                fd2 = s.Solve(net, f,delegate(IArc a, IFlowPattern _fd, double load)
                {                    
                    double sum = 0;
                    foreach (SideConstraint sc in SCList)
                    {
                        sum += Sigema * sc[a](_fd);
                    }
                    return a.PerformaceFunction(load) + sum;
                }, fd) as FlowPattern;

                //触发单次迭代结束事件
                //RasieALSIterationFinishedEvent(iter, CalculateTotalPenalty(SCList,fd2,lm_u), CalculateTotalPenaltyEvolution(SCList,lm_s, lm_u, fd, fd2));

                //foreach (SideConstraint sc in SCList)
                //{
                //    if (lm_u[sc] > 0)
                //    {
                //        lm_s[sc] = sc.LMF(fd2);
                //    }
                //    else
                //    {
                //        lm_s[sc] = Math.Max(0,Math.Pow(sc.LMF(fd2),2));
                //    }
                //}

                //④搜索步长 这里处理的是线性的边约束，如果是非线性的需要修改以下代码。
                //double Alpha = 0.5;
                //foreach (SideConstraint sc in SCList)
                //{
                //    if (lm_s[sc] < 0)
                //    {
                //        Alpha = Math.Min(Alpha, lm_u[sc] / (-1 * lm_s[sc]));
                //    }

                //}
                //double x = Alpha;

                //④搜索步长
                //double x = LS.Search(delegate (double alpha)
                //{
                //    FlowAssignment temp = s.Solve(net, f, delegate (IArc a, double load)
                //    {
                //        double sum = 0;
                //        foreach (SideConstraint sc in SCList)
                //        {
                //            sum += (lm_u[sc] + alpha * lm_s[sc]) * sc[a](fd2);
                //        }
                //        return a.PerformaceFunction(load) + sum;
                //    }, fd2) as FlowAssignment;
                //    double k = 0;
                //    foreach (IArc a in net.ArcList)
                //    {
                //        k += a.PerformaceFunction(temp.getFlow(a));
                //        foreach (SideConstraint sc in SCList)
                //        {
                //            k += lm_s[sc] * alpha;
                //        }
                //    }
                //    return k;
                //}, 0, Alpha);


                //Console.WriteLine("步长为:{0}",x);
                foreach (IArc a in net.ArcList)
                {
                    Console.WriteLine("弧{0}的流量为：{1}",a.ToString(), fd2.getFlow(a));
                }
                //foreach (SideConstraint sc in SCList)
                //{
                //    lm_u[sc] = lm_u[sc] + x * lm_s[sc];
                //    Console.WriteLine(sc.ToString()+":"+lm_u[sc]);
                //}

            }
            return fd2;
        }

        //更新拉格朗日乘子
        private static void UpdateLM(List<SideConstraint> list,SCLagrangeMultiplier lm, SCLagrangeMultiplier lm2)
        {
            foreach (SideConstraint sc in list)
            {
                lm[sc] = lm2[sc];
            }
        }
        //判断每条弧罚值是否超限
        private bool IsPenaltyLarge(List<SideConstraint> sclist, FlowPattern f, SCLagrangeMultiplier lm)
        {
            foreach (SideConstraint sc in sclist)
            {
                if (calculatePenalty(sc, f,lm) > Epsilon)
                {
                    return true;
                }
            }
            return false;
        }
        //计算总体罚值
        private double CalculateTotalPenalty(List<SideConstraint> sclist, FlowPattern f, SCLagrangeMultiplier lm)
        {
            double sum = 0;
            foreach (SideConstraint sc in sclist)
            {
                    sum += calculatePenalty(sc, f,lm);
            }
            return sum;
        }
        //计算每个边约束的罚值
        private double calculatePenalty(SideConstraint sc, FlowPattern f,SCLagrangeMultiplier lm)
        {
            double k = Math.Abs(sc.LMF(f) + (Math.Max(0, -1*Sigema * sc.LMF(f) - lm[sc]) / Sigema));
            return k;
        }
        //判断是否需要增大sigema
        private bool NeedtoEnlargeStep(List<SideConstraint> sclist, SCLagrangeMultiplier lm, SCLagrangeMultiplier lm2, FlowPattern f1, FlowPattern f2)
        {
            foreach (SideConstraint sc in sclist)
            {
                if (calculatePenalty(sc, f1, lm) != 0 && calculatePenalty(sc, f2, lm2) / calculatePenalty(sc, f1, lm) > Beta)//防止除0错误，如果一条弧k-1次迭代中没有惩罚而k次有惩罚，只能说明方向不对，再增加sigama也没有用。
                {
                    return false;
                }
            }
            return true;
        }
        //计算总体罚值改进
        private double CalculateTotalPenaltyEvolution(List<SideConstraint> sclist, SCLagrangeMultiplier lm, SCLagrangeMultiplier lm2, FlowPattern f1, FlowPattern f2)
        {
            if (CalculateTotalPenalty(sclist, f1,lm) != 0)
            {
                return CalculateTotalPenalty(sclist,f2, lm2) / CalculateTotalPenalty(sclist,f1,lm);
            }
            else
            {
                return -1;
            }
        }
    }
}
